rm(list=ls(all=TRUE))

#Set working directory to the folder with Replication Files
setwd("~/Dropbox/Nominations/Final Submission and Replication Files")

## uncomment to install packages (including dependencies) from source code OR find code below to install from CRAN
# install.packages("packages/readstata13_0.9.2.tar.gz", repo = NULL, type="source")
# install.packages("packages/readxl-1.1.0.tar.gz", repo = NULL, type="source")
# install.packages("packages/ggrepel_0.8.0.tar.gz", repo = NULL, type="source")
# install.packages("packages/gridExtra_2.3.tar.gz", repo = NULL, type="source")
# install.packages("packages/dplyr-0.7.8.tar.gz", repo = NULL, type="source")
# install.packages("packages/tidyr-0.8.2.tar.gz", repo = NULL, type="source")

#Load installed packages and/or uncomment to install from CRAN 
# install.packages("readstata13", dependencies=TRUE)
library(readstata13)
# install.packages("readxl", dependencies=TRUE)
library(readxl)
# install.packages("ggrepel", dependencies=TRUE)
library(ggrepel)
# install.packages("gridExtra", dependencies=TRUE)
library(gridExtra)
# install.packages("dplyr", dependencies=TRUE)
library(dplyr)
# install.packages("tidyr", dependencies=TRUE)
library(tidyr)

#READ IN DATA FILES
#nominee-level file 
nom.data <-  read.dta13("master_nominee_data_post_1930.dta", generate.factors=FALSE, convert.underscore = TRUE)
#short-list data
slist.data <- read.dta13("full_short_list_data_for_characteristics.dta", convert.underscore=TRUE, generate.factors=FALSE)
#dataset with factor score of presidential interest
pca.data <- read.dta13("pca_pres_interest_score.dta", convert.underscore=TRUE)
#R workspace with several summary stats of Courts of Appeals judges over time
load("COA_characteristics_workspace.Rdata")
#dataset of presidential party control (for graphing purposes)
pres.dates <- read_excel('presidential_start_stop_dates.xlsx')
#platform data
platform <- read.dta13("combined_platform_data.dta" , convert.underscore=FALSE)
#presidential speeches data
speeches <- read.dta13("updated_presidential_mentions_for_time_series_graph.dta")
#Supreme Court Central Tendency Data
sc.central.tendency <-   read.dta13("combined_central_tendency_stats_1930_2018.dta",nonint.factors= FALSE, convert.underscore=TRUE, generate.factors=FALSE)
#Supreme Court Medians at Time of Vacancy
sc.medians.at.vacancy <- read.dta13("SC_means_medians_at_time_of_vacancy_Over_Time.dta",  nonint.factors=TRUE, convert.underscore=TRUE, generate.factors=T)
#Pres-Senate Ideal Points
pres.senate <-  read.dta13("pres_senate_NOMINATE_data.dta", convert.underscore=TRUE)
#party control of Senate
party.control <- read.dta13 ("party_control_updated_1861_2018.dta")
#Post-tobit estimates from Stata
post.estimates <- read.dta13("tobit_postestimates_for_R.dta", convert.underscore = TRUE)
#IDEOLOGY SIMULATIONS
ideo.sims <-  read.dta13("ideology_sims.dta", convert.underscore = TRUE)
#DIVERSITY SIMS
div.sims <-  read.dta13("diversity_sims_from_logit.dta", convert.underscore = TRUE)

############################################################################################################
####################################
#CODE/shortcuts for graphs
########################################
cat_labeller <- function(variable,value){ return(cat_names[value])}

title.size <- 20
nom.label.size <- 10
x.axis.tick.size <- 14
y.axis.tick.size <- 14
x.axis.label.size <- 18
y.axis.label.size <- 18
strip.text.size <- 14
line.width <- 1.5
axis.text <- 1.2
label.text <- 1.2
y.offset <- .5
text.size <- 1.4
text.size2 <- 1.2
arrow.width <- 1.5
width <-2
text.size <- .7
y.axis.max <- 70
grep.size <- 4
selected.color <- "dark green"
not.color <- "purple"
gop.color <- "red"
dem.color <- "blue"

#############################################################################
#code for shading of democratic presidents in some graphs
year <- c(1893:2018)
dempres <- data.frame( matrix(data=c(year=year, dem=rep(NA, length(year))), ncol=2, nrow=length(year) ) )
names(dempres)=c("year", "dem")
dempres$dem <- ifelse(dempres$year %in% c(1893:1897, 1913:1921,1933:1953,1961:1969, 1977:1981, 1993:2001,2009:2017) , 1, 0)
out <- 7
start <- rep(NA, out)
end <- rep(NA, out)
dempres <- data.frame( matrix(data=c(start, end), nrow=out, ncol=2) )
names(dempres)=c("start", "end")
dempres$start <- c(1893, 1913,1933,1961,1977, 1993,2009)
dempres$end <- c(1897, 1921,1953,1969,1981, 2001,2017)
#pull out start year here
dempres <- subset(dempres, start>=1933)
##################################################################
############################################################################################################

############################################################################################################
#MAKE FIGURE 1
############################################################################################################
#merge PCA data with master nom data to distingish nomination years from non-nomination years
foo <- dplyr::select(nom.data,year.announced)
foo$nom.year <- 1
pca.data$year.announced <- pca.data$year #to enable merging
expanded.for.factor <- left_join(pca.data , foo)
expanded.for.factor$nom.year[is.na(expanded.for.factor$nom.year)] <- 0

pdf("Figure_1.pdf", height=6, width=11)

ggplot(data=expanded.for.factor) +
  geom_point(data=subset(expanded.for.factor, nom.year==1), aes(x=year.announced,y=pres.interest.pca), shape=19, size=3, color="red") +#years with nominations
  geom_point(data=subset(expanded.for.factor, nom.year==0), aes(x=year.announced,y=pres.interest.pca), shape=1, size=3) + #years without nominatons
  geom_line(aes(x=year.announced,y=pres.interest.pca)) + 
  geom_smooth(aes(x=year.announced,y=pres.interest.pca),  linetype="dashed", span=.5) +
  ylab("Presidential interest in Supreme Court") +
  scale_x_discrete(limits = seq(1928,2016,4)) +  theme_bw() +
  theme(        plot.title = element_text(size =title.size, face = "bold", hjust=.5),
                axis.text.x = element_text(size=x.axis.tick.size),
                axis.title.x=element_blank(),
                axis.text.y = element_text(size=y.axis.tick.size),
                axis.title.y = element_text(size=y.axis.label.size),
                panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_rect(data=dempres, aes(xmin=start-.5, xmax=end-.5, ymin=-Inf, ymax=+Inf), fill='gray20', alpha=0.2)

dev.off()

############################################################################################################
#FIGURE 2: changing nature of nominees professional background
#do separately for nominees and not
background.notnom <- filter(slist.data, nominated ==0)
background.notnom <- dplyr::select(background.notnom , date.announced, politician.cand, exec.binary.cand,ebl.binary.cand,fed.judge.cand,lawprof.binary.cand, top.law.cand, supertech.cand)
background.long.notnom <- gather(background.notnom ,category,value,  -date.announced)
background.long.notnom$nominated <- 0

background.yesnom <- filter(slist.data, nominated ==1)#just use up to Bush2 nominees for now
background.yesnom <- dplyr::select(background.yesnom , date.announced,politician.cand, exec.binary.cand,ebl.binary.cand,fed.judge.cand,lawprof.binary.cand,top.law.cand,  supertech.cand)
background.long.yesnom <- gather(background.yesnom ,category,value,  -date.announced)
background.long.yesnom$nominated <- 1
background.long <- bind_rows(background.long.notnom, background.long.yesnom)
background.long$category <- factor(background.long$category )
new.levels <- c("politician.cand","exec.binary.cand", "fed.judge.cand", "ebl.binary.cand",  "lawprof.binary.cand","top.law.cand",  "supertech.cand")
background.long$category <- factor(background.long$category, levels=new.levels )

cat_names <- list(
  'politician.cand'="Politician",
  'exec.binary.cand' = "Executive administrator",
  'fed.judge.cand'="Federal judge",
  'ebl.binary.cand'="Exec. branch lawyer",
  'lawprof.binary.cand' = "Law professor",
  'top.law.cand' = "Top law school",
  'supertech.cand'="Super Tech")

#create dummy data frame for labeling one facet
dat_text1 <- data.frame(
  label = c("","Nominees","","","", "",""), category   = unique(background.long$category),
  x     = as.Date(c(NA, "1970/1/1",NA,NA,NA,NA,NA),origin = "1970-01-01"),
  y    = c(NA, .7,NA,NA,NA,NA,NA))
dat_text2 <- data.frame(
  label = c("","Short listers","","","", "",""), category   = unique(background.long$category),
  x     = as.Date(c(NA, "1970/1/1",NA,NA,NA,NA,NA),origin = "1970-01-01"),
  y    = c(NA, .15,NA,NA,NA,NA,NA))

pdf("Figure_2.pdf", height=9, width=14)

fspan <- 1
ggplot(background.long, aes(date.announced, value)) +
  geom_smooth(data=subset(background.long, nominated==0), aes(date.announced, value), color = not.color, span=fspan, linetype = "dashed",se=FALSE)   +
  geom_smooth(data=subset(background.long, nominated==1), aes(date.announced, value), color = selected.color, span=fspan, se=FALSE) +
  geom_rug(data=subset(background.long, value == 1), aes(x = date.announced), inherit.aes = F, sides="t") + #do rug separately
  geom_rug(data=subset(background.long, value == 0), aes(x = date.announced), inherit.aes = F, sides="b") +#for top and botto
  #facet_wrap(~ category)
  facet_wrap(~ category,  labeller=cat_labeller, nrow=2,ncol=4) +
  coord_cartesian(ylim = c(0,1)) +
  geom_text(data= dat_text1, mapping = aes(x = x, y = y, label = label), color = selected.color, size=6)+
  geom_text(data = dat_text2, mapping = aes(x = x, y = y, label = label), color = not.color, size=6) +
  ylab("Probability of each type at a given point in time") + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),   
        axis.title.y = element_text(size=y.axis.label.size), 
        strip.text.x   = element_text(size = strip.text.size, colour = "dark red", angle = 0),
        plot.margin = unit(c(.5,1,.5,.5), "cm"),
        panel.spacing.x=unit(1, "lines"), panel.spacing.y=unit(.2, "lines")) 

dev.off()

########
#FIGURE 3: 
#add polygon in ggplot GOP for dem presidents
pres.dates <- filter(pres.dates,!is.na(date_start))
#kuldge fdr --> truman &  kennedy johnson
pres.dates$date_end[pres.dates$president=="roosevelt_f"] <- pres.dates$date_end[pres.dates$president=="truman"]
pres.dates <- filter(pres.dates,president!="truman")
pres.dates$date_end[pres.dates$president=="kennedy"] <- pres.dates$date_end[pres.dates$president=="johnson"]
pres.dates <- filter(pres.dates,president!="johnson")
year.start <- as.numeric(substring(pres.dates$date_start, 1,4))
year.left <- as.numeric(substring(pres.dates$date_end, 1,4))

dem.foo <- pres.dates$dem_pres
dem.data <- data.frame(year.start, year.left, dem.pres=dem.foo)
dem.data <- filter(dem.data,dem.pres==1)
dem.data <- subset(dem.data,year.start>1930)#otherwise screws up x-axis

levels <- subset(by.party.year.stats.long, grepl(".level", by.party.year.stats.long$category))
levels$Party <- factor(with(levels, ifelse(dem==1, "Democrats", "Republicans")))
levels$category <- droplevels(factor(levels$category ))
levels <- filter(levels, category!="coa.reliability.scale.levels")
levels$category <- droplevels(levels$category)
new.levels <- c("coa.exec.branch.exp.levels","coa.law.prof.levels", 
                "coa.top14.law.school.levels","coa.reliability.indicator.levels", "coa.females.levels", "coa.minorities.levels", "coa.non.white.males.levels")
levels$category <- factor(levels$category , levels = new.levels)

cat_names <- list(
  'coa.exec.branch.exp.levels'="Executive branch lawyer",
  'coa.law.prof.levels'="Law professor",
  'coa.top14.law.school.levels'="Top law school",
  'coa.reliability.indicator.levels'="EB or Law Prof. or Top school",
  'coa.females.levels'="Females",
  'coa.minorities.levels'="Minorities",
  'coa.levels.non.white.males.levels'="Female or minority")

pdf("Figure_3.pdf", height=12, width=10)

levels <- filter(levels,year>=1930)#just use post-1930

ggplot(levels, aes(x=year, y=value, color=Party, linetype=Party)) +
  geom_line(aes(group=dem)) +  
  facet_wrap(~ category,  ncol=2, labeller=cat_labeller,dir="v", scale="fixed") +
  scale_linetype_manual(values=c(1,2)) +  scale_color_manual(values=c( "blue", "red")) +
  #ggtitle("The characteristic distribution of Courts of Appeals judges, 1930-2018")  + 
  ylab("Number of each type, among all judges of a given party") + 
  theme_bw() + 
  theme( legend.position = c(0.75, 0.12), legend.text=element_text(size=18), legend.title=element_blank(),legend.direction = "horizontal", #LEGEND PARS
         panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
         plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
         axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
         axis.text.y = element_text(size=y.axis.tick.size),   
         axis.title.y = element_text(size=y.axis.label.size), 
         strip.text.x   = element_text(size = strip.text.size, colour = "dark green", angle = 0)) +
  coord_cartesian(xlim = c(1930,2018)) +
  scale_x_continuous(breaks = seq(1930,2015,15)) + 
  scale_y_continuous(breaks = seq(0,50,10)) + coord_cartesian(ylim = c(0,max(levels$value))) +
  annotate("rect", xmin=dem.data$year.start, xmax=dem.data$year.left-1, ymin=0, ymax=Inf, fill="gray80", alpha=0.5) 

dev.off()

##################################################################
#FIGURE 4: DIVERSITY OVER TIME
##################################################################
background.notnom <- filter(slist.data, nominated ==0)
background.notnom <- dplyr::select(background.notnom , date.announced,nonwhite.cand, female.cand, white.male.cand)
background.long.notnom <- gather(background.notnom ,category,value,  -date.announced)
background.long.notnom$nominated <- 0

background.yesnom <- filter(slist.data, nominated ==1)#just use up to Bush2 nominees for now
background.yesnom <- dplyr::select(background.yesnom ,  date.announced,nonwhite.cand, female.cand, white.male.cand)
background.long.yesnom <- gather(background.yesnom ,category,value,  -date.announced)
background.long.yesnom$nominated <- 1
background.long <- bind_rows(background.long.notnom, background.long.yesnom)
background.long$category <- factor(background.long$category )
new.levels <- c("white.male.cand","female.cand", "nonwhite.cand")
background.long$category <- factor(background.long$category, levels=new.levels )

cat_names <- list(
  'white.male.cand'="White males",
  'female.cand'="Females",
  'nonwhite.cand'="Minorities")

dat_text1 <- data.frame(
  label = c("Nominees","",""), category   = levels(background.long$category),
  x     = as.Date(c("2005/1/1",NA,NA),origin = "1970-01-01"),  y    = c(.75,NA,NA))
dat_text2 <- data.frame(
  label = c("Short listers","",""), category   = levels(background.long$category),
  x     = as.Date(c("1993/1/1",NA,NA),origin = "1970-01-01"),  y    = c(.55,NA,NA))

pdf("Figure_4.pdf", height=5, width=12)

fspan <- 1
ggplot(background.long, aes(date.announced, value)) +
  geom_smooth(data=subset(background.long, nominated==0), aes(date.announced, value), color = not.color, span=fspan, linetype = "dashed", se=FALSE)   +
  geom_smooth(data=subset(background.long, nominated==1), aes(date.announced, value), color = selected.color, span=fspan,se=FALSE) +
  geom_rug(data=subset(background.long, value == 1), aes(x = date.announced), inherit.aes = F, sides="t") + #do rug separately
  geom_rug(data=subset(background.long, value == 0), aes(x = date.announced), inherit.aes = F, sides="b") +#for top and botto
  facet_wrap(~ category,labeller=cat_labeller) +
  coord_cartesian(ylim = c(0,1)) +
  geom_text(data= dat_text1, mapping = aes(x = x, y = y, label = label), color = selected.color, size=4)+
  geom_text(data = dat_text2, mapping = aes(x = x, y = y, label = label), color = not.color, size=4,angle=-45) +
  ylab("Probability of each type") + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),   
        axis.title.y = element_text(size=y.axis.label.size), 
        strip.text.x   = element_text(size = strip.text.size, colour = "dark red", angle = 0),
        plot.margin = unit(c(.5,1,.5,.5), "cm"),
        panel.spacing.x=unit(1, "lines"), panel.spacing.y=unit(.2, "lines")) 

dev.off()

##################################################################
#FIGURE 5: Presidential Interest in Diversity
##################################################################
pdf("Figure_5.pdf", height=5, width=11)
ggplot(platform, aes(year, diverse_overall_index,color=party)) +
  geom_line(aes(linetype=party)) +  geom_point(aes(shape=party)) +
  scale_color_manual(values=c("blue","red")) +  scale_shape_manual(values=c(19,1)) +
  ylab("Index of diversity statements")  +
  scale_x_discrete(limits = seq(1928,2016,4)) + 
  scale_y_discrete(limits = seq(0,10,1), expand = c(0, 0)) + 
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=11),
        axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),
        axis.title.y = element_text(size=y.axis.label.size),
        legend.position = "none") +
  annotate("text", x = 1996, y =1.7, label = "Democrats",cex=5.2, col="blue") +
  annotate("text", x = 1976, y =.5, label = "Republicans",cex=5.2, col="red", angle=-75)
dev.off()

##################################################################
#FIGURE 6: Probability of non-ideological candidate
##################################################################
#make graph of predicted probability of non-censored observation for each of 3 cats
foo <- dplyr::select(post.estimates, date.announced, ideo.policy.prob.from.sl.model, rel.prob.from.nom.model, div.prob.from.nom.model, nominated,nom.name)
#keep just nominees
foo <- filter(foo, nominated==1)
foo$nominated <- NULL
foo.long <- gather(foo, category,value,  -date.announced, -nom.name)#reshape
new.levels <- c("ideo.policy.prob.from.sl.model","rel.prob.from.nom.model", "div.prob.from.nom.model")
foo.long$category <- factor(foo.long$category, levels=new.levels )
cat_names <- list('ideo.policy.prob.from.sl.model'="Ideology",'rel.prob.from.nom.model'="Reliability",'div.prob.from.nom.model'="Diversity")

pdf("Figure_6.pdf", height=5, width=12)

#for DEM shading
year <- c(1893:2018)
dempres <- data.frame( matrix(data=c(year=year, dem=rep(NA, length(year))), ncol=2, nrow=length(year) ) )
names(dempres)=c("year", "dem")
dempres$dem <- ifelse(dempres$year %in% c(1893:1897, 1913:1921,1933:1953,1961:1969, 1977:1981, 1993:2001,2009:2017) , 1, 0)
out <- 7
start <- rep(NA, out)
end <- rep(NA, out)
dempres <- data.frame( matrix(data=c(start, end), nrow=out, ncol=2) )
names(dempres)=c("start", "end")
dempres$start <- c(1893, 1913,1933,1961,1977, 1993,2009)
dempres$end <- c(1897, 1921,1953,1969,1981, 2001,2017)
#pull out start year here
dempres <- subset(dempres, start>=1933)
#convert to dates
dempres$start <- as.Date(ISOdate(dempres$start, 1, 1))
dempres$end <- as.Date(ISOdate(dempres$end, 1, 1))

ggplot(foo) +
  geom_point(data=foo, aes(date.announced, ideo.policy.prob.from.sl.model)) +
  geom_smooth(data=foo, aes(date.announced, ideo.policy.prob.from.sl.model), span=.35)  +
  #ggtitle("The probability of a non-default ideological candidate, 1930-2018")  +
  ylab("Probability of a non-default candidate") +
  theme(      plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
              axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
              axis.text.y = element_text(size=y.axis.tick.size),   
              axis.title.y = element_text(size=y.axis.label.size), 
              strip.text.x   = element_text(size = strip.text.size, colour = "dark red", angle = 0),
              plot.margin = unit(c(.5,1,.5,.5), "cm")) +
  coord_cartesian(ylim = c(.4,1)) +
  scale_y_continuous(breaks = seq(.4,1,.1))
  geom_rect(data=dempres, aes(xmin=start, xmax=end, ymin=-Inf, ymax=+Inf), fill='gray20', alpha=0.2)

dev.off()


#########################################################################################################################################################
#COUNTERFACTUAL CALCULATIONS FOR TEXT IN SECTION "Using the models to interpret changes in selection politics"
#note that every sentence in the paper that describes a calculation/statistic
#is reprinted right above the relevant calculation in the code
#########################################################################################################################################################
#########################################################################################################################################################


###################################################
#The Cost-Benefit Ratio and Nominee Ideology (Truman/Obama)
####################################
truman <- filter(post.estimates, president=="truman")

truman.actual <- truman  %>%   
  group_by(president) %>%   
  summarise(
    pi.x.adj = mean(pres.interest.pca),
    abs.president.adj = mean(abs.president),
    x.bar.adj = mean(x.bar),
    x.0 = mean(x.0),
    w.over.pi.x.adj = mean( (w.x.copar / (pi.x.adj + .01)))
  )
truman.actual
#what his mean cost-benefit ratio?
truman.actual$w.over.pi.x.adj

#now get actual predictions
truman.actual.xstar.preds <- ideo.sims$model.b.cons +
  truman.actual$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  truman.actual$abs.president.adj*ideo.sims$model.b.abs.president +
  truman.actual$x.bar.adj*ideo.sims$model.b.x.bar  +
  truman.actual$x.0*ideo.sims$model.b.x.0  
#account for truncation
truman.actual.xstar.preds <- ifelse(truman.actual.xstar.preds<0,0, truman.actual.xstar.preds)
truman.actual.xhat.preds <-   truman.actual$x.0 - truman.actual.xstar.preds

# Based on Truman's actual cost-benefit ratio, we estimate that
# for his candidates, he selected .07 units ($x_i^{\ast}$) of additional ideology in
# NOMINATE space (95\% confidence interval of [.00, .15]); combined with $x_i^0$, his
# average nominee ideology $\hat{x_i}$ was -.21 [-.29, -.14].\footnote{Because Truman is a
# Democrat, we subtract $x_i^{\ast}$ from $x_i^0$ (which takes on the average value of -.14
# for Truman) to create $\hat{x_i}$. }
truman.actual$x.0
quantile(truman.actual.xstar.preds, c(.025,.5,.975))
quantile(truman.actual.xhat.preds, c(.025,.5,.975))

#now give him Obama's ratio
truman.obama <- truman.actual
truman.obama$w.over.pi.x.adj <- with(post.estimates, mean(w.over.pi.x.copar[president =="obama"]))
truman.obama$w.over.pi.x.adj

truman.obama.xstar.preds <- ideo.sims$model.b.cons +
  truman.obama$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  truman.obama$abs.president.adj*ideo.sims$model.b.abs.president +
  truman.obama$x.bar.adj*ideo.sims$model.b.x.bar  +
  truman.obama$x.0*ideo.sims$model.b.x.0  
truman.obama.xstar.preds <- ifelse(truman.obama.xstar.preds<0,0, truman.obama.xstar.preds)
truman.obama.xhat.preds <-   truman.obama$x.0 - truman.obama.xstar.preds

#"Given Obama's cost-benefit ratio, we estimate that
# Truman would have chosen .17 [.10, .24] units of additional ideology, for a total
# ideology of -.31 [-.37, -.24].""

quantile(truman.obama.xstar.preds, c(.025,.5,.975))
quantile(truman.obama.xhat.preds, c(.025,.5,.975))

# The .09 difference between the historical high cost
# scenario and the counterfactual low-cost scenario is statistically significant [.03, .16].
quantile(truman.obama.xhat.preds-truman.actual.xhat.preds, c(.025,.5,.975))
quantile(truman.obama.xstar.preds-truman.actual.xstar.preds, c(.025,.5,.975))

# This difference is not enormous, though it amounts to one-fourth of a standard deviation
# in candidate ideology in our full sample of candidates.
quantile(truman.obama.xstar.preds-truman.actual.xstar.preds, c(.5))/sd(post.estimates$cand.dwnom.new)

################################################
#The Surprising Impact of Success: Moving the Court toward the President
################################################
with(post.estimates, mean(x.bar[year.announced == 1937]))
with(post.estimates, mean(x.bar[year.announced == 1943]))

#Use FDR: take predicted ideology in first year given actual distance, then 
#do counter-factual of last year
foo.1937 <- filter(post.estimates, nominee=="black" & nominated==1) 
#use actual values [use redundant code to keep pred code the same]
foo.1937$pi.x.adj <- with(foo.1937, mean(pres.interest.pca))
foo.1937$abs.president.adj <- with(foo.1937, mean(abs.president))
foo.1937$x.bar.adj <- with(foo.1937, mean(x.bar))
foo.1937$w.over.pi.x.adj <- with(foo.1937, (w.x.copar / (pi.x.adj + .01)))

#now get 1937 predictions
xstar.preds.1937 <- ideo.sims$model.b.cons +
  foo.1937$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  foo.1937$abs.president.adj*ideo.sims$model.b.abs.president +
  foo.1937$x.bar.adj*ideo.sims$model.b.x.bar  +
  foo.1937$x.0*ideo.sims$model.b.x.0  
xstar.preds.1937
xhat.preds.1937 <-   foo.1937$x.0  - xstar.preds.1937

# "In this scenario, the model predicts that FDR chose
# .24 [.15,.33] units of additional ideology, for a total level of -.17 [-.26, -.08]
# ($x_i^0$ was .06).
quantile(xstar.preds.1937, c(.025,.5,.975))
quantile(xhat.preds.1937, c(.025,.5,.975))

#now adjust with 1943
foo.1943 <- foo.1937
foo.1943$x.bar.adj <-  with(post.estimates, mean(x.bar[year.announced == 1943]))

xstar.preds.1943 <- ideo.sims$model.b.cons +
  foo.1943$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  foo.1943$abs.president.adj*ideo.sims$model.b.abs.president +
  foo.1943$x.bar.adj*ideo.sims$model.b.x.bar  +
  foo.1943$x.0*ideo.sims$model.b.x.0  
xstar.preds.1943
xhat.preds.1943 <-   foo.1943$x.0  - xstar.preds.1943
mean(xhat.preds.1943)
# 
# "If we substitute this level of $\bar{x}$, the predicted $x_i^{\ast}$
# and $\hat{x_i}$ fall to .18 [.08, .29] and -.12 [-.22, -.01], respectively."
quantile(xstar.preds.1943, c(.025,.5,.975))
quantile(xhat.preds.1943, c(.025,.5,.975))

#take differences
foo.1943$x.0-foo.1937$x.0

# "The .06
# difference between the actual and counter-factual is statistically different from zero
# [.02, .09]."
quantile(xhat.preds.1943-xhat.preds.1937, c(.025,.5,.975))
quantile(xstar.preds.1943-xstar.preds.1937, c(.025,.5,.975))

#how much does x-bar change?
# "By 1943, when he faced the vacancy that eventually led to the
# nomination of Rutledge, FDR had made eight appointments, shifting the mean of the Court
# from -.03 to -.20."
with(post.estimates, mean(noms.sc.mean.NSP[year.announced==1937]))
with(post.estimates, mean(noms.sc.mean.NSP[year.announced==1943]))

#does mean change before and after 1940?
foo <- filter(post.estimates, nominated==1 & president=="fdr")
# "The mean ideology of his first five picks was -.26; the mean
# of his last four nominees was -.17, noticeably more conservative."
mean(foo$cand.dwnom.new[foo$year.announced<=1940])
mean(foo$cand.dwnom.new[foo$year.announced>1940])

##################################################################
#The Importance of Seeding the Lower Courts
############################################################################
#Nixon
foo.1971 <- filter(post.estimates, nominee=="rehnquist1" & nominated==1) 
#use actual values [use redundant code to keep pred code the same]
foo.1971$pi.y.adj <- .33 #adjust interest in diversity b/c actual measure is zero; but we have private info
foo.1971$y.bar.adj <- with(foo.1971, mean(y.bar))
foo.1971$w.y.adj <- with(foo.1971, 1/ (.01+coa.nwm.copartisans))
foo.1971$w.over.pi.y.adj <- with(foo.1971, pi.y.adj / w.y.adj)

#what is probably of a NONWhite male
foo.1971.preds <- arm::invlogit( div.sims$y.hat.b.cons +
                                   foo.1971$y.bar.adj*div.sims$y.hat.b.y.bar + foo.1971$w.over.pi.y.adj*div.sims$y.hat.b.w.over.pi.y)
summary((foo.1971.preds))
# Based on the benefit-cost ratio in 1971, we estimate that Nixon had
# .07 [.03, .19] probability of nominating a non-white-male candidate.
quantile(foo.1971.preds,  c(.025,.5,.975))

#now--what if he faced same costs as Bush 2?
foo <- unique(with(post.estimates, coa.nwm.copartisans[president=="bush2"]))
foo.1971$w.y.adj.counter.bush <-  1/ (.01+foo)
foo.1971$w.over.pi.y.adj.counter.bush <- with(foo.1971, pi.y.adj /foo.1971$w.y.adj.counter.bush)

foo.1971$w.over.pi.y.adj
foo.1971$w.over.pi.y.adj.counter.bush

foo.counter.bush.preds <- arm::invlogit( div.sims$y.hat.b.cons +
                                           foo.1971$y.bar.adj*div.sims$y.hat.b.y.bar + 
                                           foo.1971$w.over.pi.y.adj.counter.bush*div.sims$y.hat.b.w.over.pi.y)
summary((foo.counter.bush.preds))
# "What if Nixon had  faced the reduced costs that George W. Bush faced in 2005 (the year he made his
# selections to the Court)? We estimate the the probability of a diverse Nixon appointee would have increased to .22 [0.06, .56]."
quantile(foo.counter.bush.preds,  c(.025,.5,.975))

#look at differences
hist (foo.counter.bush.preds-foo.1971.preds)
# "This difference of .14 is substantively meaningful and has a 95\% confidence interval of [-.01,.46]."
quantile(foo.counter.bush.preds-foo.1971.preds, c(.025,.5,.975))
mean(foo.counter.bush.preds>foo.1971.preds)

#now--what if he faced same costs as OBAMA?
foo <- unique(with(post.estimates, coa.nwm.copartisans[nominee=="sotomayor"]))
foo.1971$w.y.adj.counter.obama <-  1/ (.01+foo)
foo.1971$w.over.pi.y.adj.counter.obama <- with(foo.1971, pi.y.adj /foo.1971$w.y.adj.counter.obama)
foo.1971$w.over.pi.y.adj
foo.1971$w.over.pi.y.adj.counter.obama

foo.counter.obama.preds <- arm::invlogit( div.sims$y.hat.b.cons +
                                            foo.1971$y.bar.adj*div.sims$y.hat.b.y.bar + 
                                            foo.1971$w.over.pi.y.adj.counter.obama*div.sims$y.hat.b.w.over.pi.y)
summary((foo.counter.obama.preds))
quantile(foo.counter.obama.preds,  c(.025,.5,.975))

#look at differences
hist (foo.counter.obama.preds-foo.1971.preds)
quantile(foo.counter.obama.preds-foo.1971.preds, c(.025,.5,.975))
mean(foo.counter.obama.preds>foo.1971.preds)

# #################################################################
# #now: what if bush had faced same costs as obama [Not used in paper]
# ####################################################
# foo.2005 <- filter(post.estimates, nominee=="miers" & nominated==1) 
# #use actual values [use redundant code to keep pred code the same]
# foo.2005$pi.y.adj <- .33 #kludging some interest in diversity b/c actual measure is zero; but we have private info
# foo.2005$y.bar.adj <- with(foo.2005, mean(y.bar))
# foo.2005$w.y.adj <- with(foo.2005, 1/ (.01+coa.nwm.copartisans))
# foo.2005$w.over.pi.y.adj <- with(foo.2005, pi.y.adj / w.y.adj)
# 
# #what is probably of a NONWhite male
# foo.2005.preds <- arm::invlogit( div.sims$y.hat.b.cons +
#                                    foo.2005$y.bar.adj*div.sims$y.hat.b.y.bar + foo.2005$w.over.pi.y.adj*div.sims$y.hat.b.w.over.pi.y)
# summary((foo.2005.preds))
# quantile(foo.2005.preds,  c(.025,.5,.975))
# 
# 
# #now--what if he faced same costs as Bush 2?
# foo <- with(post.estimates, coa.nwm.copartisans[president=="bush2"])
# foo.2005$w.y.adj.counter <-  unique(1/ (.01+foo))
# foo.2005$w.over.pi.y.adj.counter <- with(foo.2005, pi.y.adj /foo.2005$w.y.adj.counter)
# 
# foo.2005$w.over.pi.y.adj
# foo.2005$w.over.pi.y.adj.counter
# 
# foo.counter.preds <- arm::invlogit( div.sims$y.hat.b.cons +
#                                       foo.2005$y.bar.adj*div.sims$y.hat.b.y.bar + 
#                                       foo.2005$w.over.pi.y.adj.counter*div.sims$y.hat.b.w.over.pi.y)
# summary((foo.counter.preds))
# quantile(foo.counter.preds,  c(.025,.5,.975))


########################################################################################################
#APPENDIX: The importance of the lower court pipeline for ideology
#################################################### ####################################################
#1--look at effect of "seeding" --compare Bush/Reagan early on in their term to late, using CO-par measure
foo.1981 <- filter(post.estimates, nominee=="oconnor" & nominated==1) 
#set covariates at mean of reagan/bush
foo.1981$pi.x.adj <- with(post.estimates, mean(pres.interest.pca[president %in% c("bush", "reagan")])) 
foo.1981$abs.president.adj <- with(post.estimates, mean(abs.president[president %in% c("bush", "reagan")])) 
foo.1981$x.bar.adj <- with(post.estimates, mean(x.bar[president %in% c("bush", "reagan")])) 
#now, let costs and w_o vary based on year  
foo.1981$w.over.pi.x.adj <- with(foo.1981, (w.x.copar / (pi.x.adj + .01)))

#now get 1981 predictions
xstar.preds.1981 <- ideo.sims$model.b.cons +
  foo.1981$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  foo.1981$abs.president.adj*ideo.sims$model.b.abs.president +
  foo.1981$x.bar.adj*ideo.sims$model.b.x.bar  +
  foo.1981$x.0*ideo.sims$model.b.x.0  
xhat.preds.1981 <-   foo.1981$x.0 +xstar.preds.1981

# "We estimate that Reagan in 1981 would have chosen .31 [.21, .40]
# additional units of ideology, for a total ideology of .17 [.08, .26]; recall the default
# level of ideology was .11."
foo.1981$x.0
quantile(xhat.preds.1981, c(.025,.5,.975))
quantile(xstar.preds.1981, c(.025,.5,.975))

#now do 1991
foo.1991 <- filter(post.estimates, nominee=="souter" & nominated==1) 
#set covariates at mean of reagan/bush
foo.1991$pi.x.adj <- with(post.estimates, mean(pres.interest.pca[president %in% c("bush", "reagan")])) 
foo.1991$abs.president.adj <- with(post.estimates, mean(c( unique(abs.president[president=="bush"]), unique(abs.president[president=="reagan"]))))
foo.1991$x.bar.adj <- with(post.estimates, mean(x.bar[president %in% c("bush", "reagan")])) 
#now, let costs and w_o vary based on year  
foo.1991$w.over.pi.x.adj <- with(foo.1991, (w.x.copar / (pi.x.adj + .01)))

#now get 1991 predictions
xstar.preds.1991 <- ideo.sims$model.b.cons +
  foo.1991$w.over.pi.x.adj*ideo.sims$model.b.w.over.pi.x +
  foo.1991$abs.president.adj*ideo.sims$model.b.abs.president +
  foo.1991$x.bar.adj*ideo.sims$model.b.x.bar  +
  foo.1991$x.0*ideo.sims$model.b.x.0  
xhat.preds.1991 <-   foo.1991$x.0 +xstar.preds.1991

# "This simulation estimates that Bush chose a
# comparable amount of additional ideology (.27, [.20,.33]) as Reagan did in 1981; however,
# because of the increase in $x_i^0$, the estimate of $\hat{x}$ for Bush rises to .4 [.33,.47]."
quantile(xhat.preds.1991, c(.025,.5,.975))
quantile(xstar.preds.1991, c(.025,.5,.975))

#take differences
foo.1991$x.0-foo.1981$x.0
par(mfrow=c(1,2))
hist (xstar.preds.1991-xstar.preds.1981)
hist (xhat.preds.1991-xhat.preds.1981)
" The .23 difference in predicted overall ideology is statistically significant [.11,
.36] and substantively large"
quantile(xhat.preds.1991-xhat.preds.1981, c(.025,.5,.975))
quantile(xstar.preds.1991-xstar.preds.1981, c(.025,.5,.975))

# #what is x-o and COA co-partisans in both years?
# foo.1981$coa.n.copartisans
# foo.1981$x.0
# foo.1991$coa.n.copartisans
# foo.1991$x.0


####################################################################################################################################
####################################################################################################################################
#APPENDIX FIGURES
####################################################################################################################################
####################################################################################################################################

##################################################################
#FIGURE A-1: Aggregate index of presidential/party interest, 1928-2016
##################################################################

pdf("Figure_A1.pdf", height=6, width=11)

ggplot(platform, aes(year, index_agg, color=party)) +
  geom_line(aes(linetype=party)) +  geom_point(aes(shape=party)) +
  scale_color_manual(values=c("blue","red")) +  scale_shape_manual(values=c(19,1))+
  ylab("Aggregate index of party/pres.\ninterest in Supreme Court")  + 
  scale_fill_manual(values=c("blue", "red")) + 
  scale_x_discrete(limits = seq(1928,2016,4)) + 
  scale_y_discrete(limits = seq(0,20,5), expand = c(0, 1)) + 
  theme(plot.margin=unit(c(.5,1,.5,.5),"lines"),
        legend.position="none",
        plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size),
        axis.title.x=element_blank(), 
        axis.text.y = element_text(size=y.axis.tick.size),
        axis.title.y = element_text(size=y.axis.label.size)) +
  annotate("text", x = 2008, y =1.5, label = "Democrats",cex=5.2, col="blue") +
  annotate("text", x = 2007, y =13, label = "Republicans",cex=5.2, col="red")

dev.off()

##################################################################
#Figure A-2: Presidential Interest in Supreme Court Policy, 1930-2018.
##################################################################

pdf ("Figure_A2.pdf", height=5, width=12)

speeches <- filter(speeches, year>=1930)
#convert to dates for shading purpposes
x <- as.Date( paste0(speeches$year,"/7/1"))
x.axis.label <- seq(1935,2015,10)
x.axis.loc <- seq.Date(as.Date("1935-07-01"), as.Date("2015-07-01"), by="10 years")

par(mfrow=c(1,1), mar=c(2,4,1,2))
axis.text <- 1.2
plot(x, speeches$presimport, type="n", axes=FALSE, xlab="", ylab="", xlim=c(as.Date("1933-01-01"), as.Date("2015-07-31")), xaxs="r")

#add shading for years in which president was democrat
shade.color <- "gray80"
polygon(x=c(as.Date("1893/03/04"), as.Date("1893/03/04"), as.Date("1897/03/04"), as.Date("1897/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,    border=F)

polygon(x=c(as.Date("1913/03/04"), as.Date("1913/03/04"), as.Date("1921/03/04"),as.Date("1921/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1933/03/04"), as.Date("1933/03/04"), as.Date("1953/01/20"), as.Date("1953/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1961/01/20"), as.Date("1961/01/20"), as.Date("1969/01/20"), as.Date("1969/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)           
polygon(x=c(as.Date("1977/01/20"), as.Date("1977/01/20"), as.Date("1981/01/20"), as.Date("1981/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)            
polygon(x=c(as.Date("1993/01/20"), as.Date("1993/01/20"), as.Date("2001/01/20"), as.Date("2001/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         
polygon(x=c(as.Date("2009/01/20"), as.Date("2009/01/20"), as.Date("2017/01/20"), as.Date("2017/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         
lines(x,  speeches$presimport, type="b", pch=19, bg="white")
axis.Date(side = 1, at=x.axis.loc, label=x.axis.label,mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, at = seq(0,14,2), mgp=c(2,.5,0), cex.axis=axis.text, las=1)
box()
mtext("Presidential speech index",2, cex=1.2, srt=90, line=2)
#add labels for each president

pres.dates <- subset(pres.dates, date_start > "1923-08-02")#keep only from hooever
pres.dates$date_start <- as.Date(pres.dates$date_start)
pres.dates$date_end <- as.Date(pres.dates$date_end)
segments(pres.dates$date_start, par()$usr[3], pres.dates$date_start,par()$usr[4],lty=2, col="pink", lwd=2)
titles <- tools::toTitleCase(pres.dates$president)
titles <- ifelse(titles == "Roosevelt_f", "FDR", 
                 ifelse( titles == "Bush", "Bush 1",  
                         ifelse( titles == "Bush2", "Bush 2",titles)))
#find midpoint of each president
mid <- pres.dates$date_start + floor((pres.dates$date_end-pres.dates$date_start)/2)
text(mid-200, par()$usr[4]-.5, titles, xpd=NA, col="blue", cex=.9, srt=90, adj=1)
box()
#add rug for each nominations
nom.data <- arrange(nom.data, date.offic.nom)
foo <- nom.data$date.offic.nom
y <- rep(par()$usr[3]+.02, length(nom.data$date.offic.nom))
text(foo+180,y, "|")

dev.off()

##################################################################
#FIGURE A3: The distribution of lower court characteristics over time.
##################################################################
rates <- subset(by.party.year.stats.long, grepl("rate", by.party.year.stats.long$category))
#pull out  coa.reliability scale becaue confusing
rates <- subset(rates, rates$category !="coa.reliability.scale.rate" & rates$category != "coa.reliability.index.rate" & rates$category != "coa.diversity.index.rate")
rates$category <- droplevels(factor(rates$category ))
n.cats <- length(unique(rates$category))
new.levels <- c("coa.exec.branch.exp.rate","coa.law.prof.rate", "coa.top14.law.school.rate","coa.reliability.indicator.rate", "coa.females.rate", "coa.minorities.rate", "coa.non.white.males.rate")
#create party indicator for graph
rates$Party <- factor(with(rates, ifelse(dem==1, "Democrats", "Republicans")))
rates$category <- factor(rates$category , levels = new.levels)
#for text
foo <- filter(rates, year==1982 & dem==0)
cat_names <- list(
  'coa.exec.branch.exp.rate'="Executive branch lawyer",
  'coa.law.prof.rate'="Law professor",
  'coa.top14.law.school.rate'="Top law school",
  'coa.reliability.indicator.rate'="Overall reliabilty indicator (EB or Law Prof. or Top school)",
  #'coa..reliability.scale.rate'="Overall reliabilty scale (EB or Law Prof. or Top school)",
  'coa.females.rate'="Females",
  'coa.minorities.rate'="Minorities",
  'coa.non.white.males.rate'="Female or minority")

pdf("Figure_A3.pdf", height=12, width=13)

rates <- filter(rates,year>=1930)
ggplot(rates, aes(x=year, y=value, color=Party, linetype=Party)) +
  geom_line(aes(group=dem)) +  
  facet_wrap(~ category,  ncol=2, labeller=cat_labeller,dir="v") +
  #facet_wrap(~ category,  ncol=2) 
  scale_linetype_manual(values=c(1,2)) +  scale_color_manual(values=c( "blue", "red")) +
  ggtitle("The characteristic distribution of Courts of Appeals judges, 1930-2018")  + 
  ylab("Proportion of each type, among all judges of a given party") + 
  theme_bw() + 
  theme( panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
         plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
         axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
         axis.text.y = element_text(size=y.axis.tick.size),   
         axis.title.y = element_text(size=y.axis.label.size), 
         strip.text.x   = element_text(size = strip.text.size, colour = "dark green", angle = 0),
         legend.position = c(0.75, 0.1), legend.text=element_text(size=18), legend.title=element_blank(), legend.direction = "horizontal") +
  scale_x_continuous(breaks = seq(1930,2015,15)) + coord_cartesian(xlim = c(1930,2018)) +
  scale_y_continuous(breaks = seq(0,1,.2)) + coord_cartesian(ylim = c(0,1)) +
  annotate("rect", xmin=dem.data$year.start, xmax=dem.data$year.left, ymin=0, ymax=Inf, fill="gray80", alpha=0.5) 

dev.off()

##################################################################
#FIGURE A4: separate plot for reliability index.
##################################################################

rel.index.party <- subset(by.party.year.stats.long, category=="coa.reliability.index.rate")
rel.index.party$party <- with(rel.index.party , ifelse(dem==0, "GOP", "DEM"))
rel.index.party$dem <- NULL
rel.index.party$category <- NULL
#add overall as well
rel.index.all  <- subset(year.level.stats, select=c(year, coa.mean.reliability.index))
rel.index.all$party <- "ALL"
rel.index.all <- rename(rel.index.all,  value="coa.mean.reliability.index")
foo <- bind_rows(rel.index.party, rel.index.all)
foo <- arrange( filter(foo, year>=1930), year, party)
#for text, what is rate ?
with(foo, value[year==1953 & party=="GOP"])
with(foo, value[year==2005 & party=="GOP"])
#now using just overall
foo <- filter(foo, party=="ALL")

pdf("Figure_A4.pdf", height=5, width=12)

ggplot(foo, aes(year, value)) + geom_line() + geom_smooth() +
  ggtitle("The mean reliability of Courts of Appeals judges")  +
  ylab("Mean reliability index") + 
  theme_bw() + 
  theme(legend.position="none",
    plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
    axis.text.x = element_text(size=x.axis.tick.size),  axis.title.x=element_blank(), 
    axis.text.y = element_text(size=y.axis.tick.size),   
    axis.title.y = element_text(size=y.axis.label.size), 
    strip.text.x   = element_text(size = strip.text.size, colour = "dark green", angle = 0)) +
  coord_cartesian(xlim = c(1930,2018)) +   scale_x_continuous(breaks = seq(1930,2015,15)) + 
  scale_y_continuous(breaks = seq(1.3,2,.1)) + coord_cartesian(ylim = c(1.3,2)) +
  annotate("rect", xmin=dem.data$year.start, xmax=dem.data$year.left-1, ymin=-Inf, ymax=Inf, fill="gray80", alpha=0.5) 

dev.off()

##################################################################
#FIGURE A5: The ideology of the president  and the existing Court’s respective ideology  over time
##################################################################

pdf("Figure_A5.pdf", height=6, width=12)

axis.text <- 1.2
title.text <- 1.2
label.text <- 1.4
line.width <- 1.5
title.size <- 20
nom.label.size <- 10
x.axis.tick.size <- 14
y.axis.tick.size <- 14
y.axis.label.size <- 18
pres.color <- "dark green"
NSP.color <- "orange"
MQ.color <- "red"

par(mar=c(2,4,3.5,1))
#first do presidnents
x <- seq(as.Date("1930/1/1"), as.Date("2018/1/1"), "years")
y.pres <- sc.central.tendency$pres.dwnom1
x.axis.label <- seq(1930,2020,5)
x.axis.loc <- seq(as.Date("1930-01-01"), as.Date("2020-07-31"), length.out=19)

plot(x, y.pres, axes=FALSE, type="n", xlab="", ylab="", xaxs="i", xlim=c(as.Date("1930-01-01"), as.Date("2020-07-31")))

#add shading for dem presidents


#add shading for years in which president was democrat
shade.color <- "gray80"
polygon(x=c(as.Date("1893/03/04"), as.Date("1893/03/04"), as.Date("1897/03/04"), as.Date("1897/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,    border=F)

polygon(x=c(as.Date("1913/03/04"), as.Date("1913/03/04"), as.Date("1921/03/04"),as.Date("1921/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1933/03/04"), as.Date("1933/03/04"), as.Date("1953/01/20"), as.Date("1953/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1961/01/20"), as.Date("1961/01/20"), as.Date("1969/01/20"), as.Date("1969/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)           
polygon(x=c(as.Date("1977/01/20"), as.Date("1977/01/20"), as.Date("1981/01/20"), as.Date("1981/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)            
polygon(x=c(as.Date("1993/01/20"), as.Date("1993/01/20"), as.Date("2001/01/20"), as.Date("2001/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         
polygon(x=c(as.Date("2009/01/20"), as.Date("2009/01/20"), as.Date("2017/01/20"), as.Date("2017/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         
axis.Date(side = 1, at=x.axis.loc, label=x.axis.label,mgp=c(2,.5,0), cex.axis=axis.text)
axis(2)
#add labels for each president
pres.dates <- subset(pres.dates, date_start > "1923-08-02")#keep only from hooever
pres.dates$date_start <- as.Date(pres.dates$date_start)
pres.dates$date_end <- as.Date(pres.dates$date_end)
segments(pres.dates$date_start, -3, pres.dates$date_start,3,lty=2, col="pink", lwd=2)
titles <- tools::toTitleCase(pres.dates$president)
titles <- ifelse(titles == "Roosevelt_f", "FDR", 
                 ifelse( titles == "Bush", "Bush 1",  
                         ifelse( titles == "Bush2", "Bush 2",titles)))
#find midpoint of each president
mid <- pres.dates$date_start + floor((pres.dates$date_end-pres.dates$date_start)/2)
text(mid, .75, titles, xpd=NA, col="blue", cex=.9, srt=45, adj=0)
box()
#add flat lines for president (ideal points don't vary)
foo <- dplyr::select(sc.central.tendency, president, pres.dwnom1)
foo <- foo[!duplicated(foo), ]
foo <- full_join(foo, pres.dates)
segments(foo$date_start,foo$pres.dwnom1, foo$date_end,foo$pres.dwnom1, col=pres.color, lwd=4)
text(as.Date("1950/1/1"), -.43,"Pres.\nideal point", col=pres.color)

#now do Court
#mean
y  <- sc.central.tendency$sc.mean.NSP[!is.na(sc.central.tendency$sc.mean.NSP)]
x  <- seq(as.Date("1930/1/1"), as.Date("2018/1/1"), "years")[!is.na(sc.central.tendency$sc.mean.NSP)]#account for missing years in NSP
points(x, y, col=NSP.color, type="l", lty=1, lwd=2)
text(as.Date("1949/1/1"), -.1, "Mean of Court (NSP)", col=NSP.color, srt=0)

#now add tick marks for nominations
sc.medians.at.vacancy$date.offic.nom
y <- rep(par()$usr[3]+.02, length(sc.medians.at.vacancy$date.offic.nom))
text(sc.medians.at.vacancy$date.offic.nom,y, "|")

dev.off()

##################################################################
#FIGURE A6: Scatterplot of NSP measure versus Nemacheck measure of nominee ideology, 1930- 2018, for Supreme Court nominees.
##################################################################

nom.data <- filter(slist.data, nominated==1 & nominee!="roberts_j_2")

pdf("Figure_A6.pdf", height=6, width=8)

nom.data <- mutate(nom.data, x =nsp1.new, y=cand.dwnom.new)
foo  <- with(nom.data, cor(x,y))
foo <- paste("Correlation = ", round(foo,2))
ggplot(nom.data, aes(x,y)) + 
  geom_text_repel(data=subset(nom.data, dem.president==0), aes(x,y, label = nom.name, fontface = "italic"), size=grep.size, color=gop.color) +
  geom_text_repel(data=subset(nom.data, dem.president==1), aes(x,y, label = nom.name), size=grep.size, color=dem.color) +
  geom_smooth(col="black") +
  xlab("NSP score (liberal to conservative)") + ylab("Nemacheck (liberal to conservative)") +
  scale_x_continuous(breaks = seq(-.75,.75,.15)) +  scale_y_continuous(breaks = seq(-.75,.75,.15)) +
  coord_cartesian(xlim = c(-75,.75)) + coord_cartesian(ylim = c(-.75,.75)) +
  theme(plot.title = element_text(size =title.size, face = "bold", hjust=.5), 
        axis.text.x = element_text(size=x.axis.tick.size), 
        axis.text.y = element_text(size=y.axis.tick.size),  
        axis.title.x = element_text(size=x.axis.label.size),
        axis.title.y = element_text(size=y.axis.label.size)) +
  annotate("text", x = .5, y = -.6, label =foo,cex=5.2)

dev.off()


##################################################################
#FIGURE A7: just number of democrats and republicans on the Courts of Appeals.
##################################################################
foo <- dplyr::select(year.level.stats,year, coa.n.dem, coa.n.gop)
foo <- filter(foo, year>=1930)
foo <- gather(foo, party, number.judges, -year)
ok <- foo$party == "coa.n.dem"

pdf("Figure_A7.pdf", height=5, width=12)

x <- as.Date( paste0(foo$year,"/7/1"))
x.axis.label <- seq(1935,2015,10)
x.axis.loc <- seq.Date(as.Date("1935-07-01"), as.Date("2015-07-01"), by="10 years")
axis.text <- 1.2

par(mfrow=c(1,1), mar=c(2,4,1,2))

plot(x, foo$number.judges, type="n", xaxs="r", axes=FALSE, xlab="", ylab="", xlim=c(as.Date("1933-01-01"), as.Date("2015-07-31")))

#add shading for years in which president was democrat
shade.color <- "gray80"
polygon(x=c(as.Date("1893/03/04"), as.Date("1893/03/04"), as.Date("1897/03/04"), as.Date("1897/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,    border=F)

polygon(x=c(as.Date("1913/03/04"), as.Date("1913/03/04"), as.Date("1921/03/04"),as.Date("1921/03/04")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1933/03/04"), as.Date("1933/03/04"), as.Date("1953/01/20"), as.Date("1953/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)        
polygon(x=c(as.Date("1961/01/20"), as.Date("1961/01/20"), as.Date("1969/01/20"), as.Date("1969/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)           
polygon(x=c(as.Date("1977/01/20"), as.Date("1977/01/20"), as.Date("1981/01/20"), as.Date("1981/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)            
polygon(x=c(as.Date("1993/01/20"), as.Date("1993/01/20"), as.Date("2001/01/20"), as.Date("2001/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         
polygon(x=c(as.Date("2009/01/20"), as.Date("2009/01/20"), as.Date("2017/01/20"), as.Date("2017/01/20")),
        y=par()$usr[c(3,4,4,3)],
        col= shade.color,   
        border=F)         

lines(x[ok], foo$number.judges[ok], col="blue")
lines(x[!ok], foo$number.judges[!ok], col="red", lty=2)
axis.Date(side = 1, at=x.axis.loc, label=x.axis.label,mgp=c(2,.5,0), cex.axis=axis.text)
axis(2, at =seq(0,110,10), mgp=c(2,.5,0), cex.axis=1.2, las=1)
mtext("Number of judges on the Courts of Appeals", 2, cex=1.5, srt=90, line=2.5)
box()
text(-6408, 60, "Democratic judges", col="blue", cex=1.4)
text(-5400, 11, "Republicans judges", col="red", cex=1.4, adj=0)
########################################################
#add labels for each president
pres.dates <- read_excel('presidential_start_stop_dates.xlsx')
pres.dates <- subset(pres.dates, date_start > "1923-08-02")#keep only from hooever
pres.dates$date_start <- as.Date(pres.dates$date_start)
pres.dates$date_end <- as.Date(pres.dates$date_end)
segments(pres.dates$date_start, par()$usr[3], pres.dates$date_start,par()$usr[4],lty=2, col="pink", lwd=2)
titles <- tools::toTitleCase(pres.dates$president)
titles <- ifelse(titles == "Roosevelt_f", "FDR", 
                 ifelse( titles == "Bush", "Bush 1",  
                         ifelse( titles == "Bush2", "Bush 2",titles)))
#find midpoint of each president
mid <- pres.dates$date_start + floor((pres.dates$date_end-pres.dates$date_start)/2)
text(mid-100, par()$usr[4]-1, titles, xpd=NA, col="blue", cex=.9, srt=90, adj=1)
box()

dev.off()

##################################################################
#FIGURE A8: The cost-benefit ratio for presidential demand of ideology
##################################################################
party.control <- dplyr::select(party.control, year, senate_dems, senate_gops)
party.control <- party.control[rep(1:nrow(party.control),each=2),] 
party.control$year <- 1860 + (1:nrow(party.control))
pca.data$year <- pca.data$year.announced

pres.senate <- filter ( dplyr::select(pres.senate, year, dem.president), year>=1928)
pres.senate <- inner_join(pres.senate,pca.data)
pres.senate <- inner_join(pres.senate, party.control)

#merge in copartisans and senate
year.data.for.ratio <- dplyr::select(year.level.stats, year, coa.n.dem, coa.n.gop)
year.data.for.ratio <- dplyr::filter(year.data.for.ratio, year>=1930)
year.data.for.ratio <- inner_join(pres.senate,year.data.for.ratio)

#now construct ratios
year.data.for.ratio <- mutate(year.data.for.ratio,
                              coa_n_copartisans = ifelse(dem.president==1, coa.n.dem,coa.n.gop),
                              pi_x = pres.interest.pca,
                              w_x_copar = 1/coa_n_copartisans,
                              w_over_pi_x_copar =  w_x_copar / (pi_x + .01))
year.data.for.ratio$w_over_pi_x_copar
#what years have zero interest
with(year.data.for.ratio, year[pi_x==0])

#exclude years with zero interest, which drives costs way up
foo <- filter(year.data.for.ratio, pi_x!=0 )

copar <- ggplot(foo,  aes(x=year, y=w_over_pi_x_copar)) +
  geom_point()+ geom_smooth(span=1)+ 
  ggtitle("The cost-benefit ratio of ideology")  +
  ylab("Ratio (based on co-partisan judges)") + 
  theme_bw() + 
  theme(        plot.title = element_text(size =title.size, face = "bold", hjust=.5),
                axis.text.x = element_text(size=x.axis.tick.size),
                axis.title.x=element_blank(),
                axis.text.y = element_text(size=y.axis.tick.size),
                axis.title.y = element_text(size=y.axis.label.size),
                panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  coord_cartesian(xlim = c(1930,2018)) +   scale_x_continuous(breaks = seq(1930,2015,15)) + 
  annotate("rect", xmin=dem.data$year.start, xmax=dem.data$year.left-1, ymin=-Inf, ymax=Inf, fill="gray80", alpha=0.5) 

#now using senate
year.data.for.ratio <- mutate(year.data.for.ratio,
                              total = senate_dems + senate_gops,
                              percent_sameparty_senators = ifelse(dem.president==1, senate_dems/total, senate_gops/total),
                              w_x_senate = 1/percent_sameparty_senators,
                              w_over_pi_x_senate =  w_x_senate / (pi_x + .01))
foo <- filter(year.data.for.ratio, pi_x!=0 )

foo <- filter(year.data.for.ratio ,w_over_pi_x_copar<2 )
sen <- ggplot(foo,  aes(x=year, y=w_over_pi_x_senate)) +
  geom_point()+ geom_smooth(span=1)+ 
  ggtitle("The cost-benefit ratio of ideology")  +
  ylab("Ratio (based on Senate )") + 
  theme_bw() + 
  theme(        plot.title = element_text(size =title.size, face = "bold", hjust=.5),
                axis.text.x = element_text(size=x.axis.tick.size),
                axis.title.x=element_blank(),
                axis.text.y = element_text(size=y.axis.tick.size),
                axis.title.y = element_text(size=y.axis.label.size),
                panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  coord_cartesian(xlim = c(1930,2018)) +   scale_x_continuous(breaks = seq(1930,2015,15)) + 
  annotate("rect", xmin=dem.data$year.start, xmax=dem.data$year.left-1, ymin=-Inf, ymax=Inf, fill="gray80", alpha=0.5) 


pdf("Figure_A8.pdf", height=5, width=15)
grid.arrange(copar , sen, nrow = 1)
dev.off()



##################################################################
#FIGURE A9: Policy reliability index, for selected candidates, 1930-2018.
##################################################################

pdf ("Figure_A9.pdf", height=8, width=10)

foo <- filter(slist.data, nominated==1) # Dotpot--just for nominees
foo <- arrange(foo, date.announced)

par(mfrow=c(1,1), mar=c(4,6,.5,1))

rel.sort <- sort(foo$reliable.cand.binary, decreasing=TRUE)
nom.sort <- foo$nom.name[order(foo$reliable.cand.binary, foo$nom.name, decreasing=TRUE)]
y <- c(1:length(nom.sort))
axis.text <- .8
plot(rel.sort,y, pch=19, axes=FALSE, type="n",xlab="", ylab="")
abline(h=y, col="gray80")
points(rel.sort,y, pch=19)
#segments(-1, y, rel.sort,y, lty=2, lwd=2)
axis(2, at=y, labels=nom.sort, las=1, mgp=c(2,.5,0), cex.axis=axis.text)
axis(1 ,mgp=c(2,.5,0), at=seq(0,4,1), cex.axis=axis.text+.2)
box()
mtext("Reliability index", 1, line=2, cex=1.2)

dev.off()


##################################################################
##################################################################
#FORMAL THEORY GRAPHS FOR APPENDIX B (THESE DO NOT INVOLVE ANY DATA)
##################################################################
##################################################################

text.size <- 1.4
y.bottom <- -.16

#illustrate nominee buy
pdf("Figure_B1A.pdf", width=7, height=2)

par(mar=c(0,0,0,0))
plot(0,0, type="n", xlim=c(-.1,1.1), ylim=c(-.2,.6),xlab="", ylab="",axes=F, main="")
arrows(0,0,1.1,0,lwd=4,font=2, code=3)
segments(.15,-.1,.15,.1,lwd=4,font=2)
text(.15,y.bottom,0, cex=text.size)
segments(.3,-.1,.3,.1,lwd=4,font=2)
text(.3,y.bottom,expression(italic(x[0])), cex=text.size)
arrows(.55,.3,.55,.01,lwd=4,font=2, lty=1, col="blue", length=.2)

text(.55,.4,expression(paste(italic(a),italic(bar(x)[i])," + ",italic(b), italic(x[i]^0))), cex=text.size)
segments(.6,-.1,.6,.1,lwd=4,font=2)
text(.6,y.bottom,expression(bar(italic(x))), cex=text.size)
segments(.8,-.1,.8,.1,lwd=4,font=2)
text(.8,y.bottom,expression(italic(p)), cex=text.size)
arrows(1,.3,1,.01,lwd=4,font=2, lty=1, col="blue", length=.2)
text(1,.45,expression(paste(italic(frac(p,b)), 
                            " - ",  italic(frac(a,b)),   italic(bar(x)[i]),
                            " - ", italic(x[i]^0))), cex=text.size)
text(.1,.47, "A) Example of\nideology choice", cex=1.5, font=3)

dev.off()

#now as a function of cost of ideology
width <- 10
axis.text <- 1.2

pdf("Figure_B1B.pdf", width=8, height=6)

par(mar=c(4,4,1,1))
plot(0,0, type="n", xlim=c(0,.25), ylim=c(0,.7),xlab="", ylab="",axes=F, main="", xaxs="i", yaxs="i")
segments(0,.7, .15,0, lwd=5, col="blue")
segments(.15,0,.25,0,  lwd=width, col="blue")
axis(1, at = seq(0,.25,.05), mgp=c(2,.5,0), las=1, cex.axis=axis.text)
axis(2, at = seq(0,.7,.1), mgp=c(2,.5,0), las=1, cex.axis=axis.text)
mtext("Cost of choosing additional ideology", 1, cex=text.size, line=2.5, font=2)
mtext("Additional ideology chosen", 2, cex=text.size, line=2.5, font=2)
box()
text(.17,.6, "B) Ideology\nagainst costs", cex=2, font=3)

dev.off()


